﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Drawing;

namespace erse253_homework
{
    public class NumericalSolutions
    {

        /// <summary>
        /// Calculates the average value of the data in the range of the vector.
        /// </summary>
        /// <param name="calc">
        /// Dataset
        /// </param>
        /// <returns>
        /// Mean
        /// </returns>
        public double Mean(double[,]  calc)
        {
            double sum = 0;
            int n = 0;
            if (calc != null)
            {
                for (int i = 0; i < calc.GetLength(0); i++)
                    for (int j = 0; j < calc.GetLength(1); j++)
                        if (!double.IsNaN(calc[i, j]))
                        {
                            sum += calc[i, j];
                            n++;
                        }

                if (n > 0)
                    return sum / (n);
                else
                    return double.NaN;


            }
            return double.NaN;
        }

        public int NumberOfPoints(double[,] calc)
        {
            int n = 0;
            int x = calc.GetLength(0);
            int y = calc.GetLength(1);


            if (calc != null)
            {
                for (int i = 0; i < x; i++)
                    for (int j = 0; j < y; j++)
                        if (!double.IsNaN(calc[i, j]))
                            n++;

            }
            return n;
        }
        public double MeanSpecial(double[,] calc)
        {
            double sum = 0;
            int n = 0;
            int x = calc.GetLength(0);
            if(x > 1) 
                x--;

            int y = calc.GetLength(1);
            if (y > 1) 
                y--;

            if (calc != null)
            {
                if (calc != null)
                    for (int i = 0; i < x; i++)
                        for (int j = 0; j < y; j++)
                            if (!double.IsNaN(calc[i, j]))
                            {
                                sum += calc[i, j];
                                n++;
                            }


                if (n > 0)
                    return sum / (n);
                else
                    return double.NaN;


            }
            return double.NaN;
        }


        /// <summary>
        /// Calculates the average value of the data in the range of the vector.
        /// </summary>
        /// <param name="calc">
        /// Dataset
        /// </param>
        /// <returns>
        /// Mean
        /// </returns>
        public double Mean(double[] calc)
        {
            double sum = 0;

            if (calc != null)
            {
                for (int i = 0; i < calc.Length; i++)
                    sum += calc[i];

                return sum / (calc.Length);
            }
            return double.NaN;
        }

        /// <summary>
        /// Calculates the median value of the data in the range of the vector.
        /// </summary>
        /// <param name="calc">
        /// Dataset
        /// </param>
        /// <returns>
        /// Median
        /// </returns>
        public double Median(double[,] calc) 
        {
            if (calc != null)
            {
                double[] array = new double[calc.Length];

                for (int i = 0; i < calc.GetLength(0); i++)
                    for (int j = 0; j < calc.GetLength(1); j++)
                        array[i + j * calc.GetLength(0)] = calc[i, j];

                Array.Sort(array);


                return ReturnMedian(array);
            }
            return double.NaN;
     
        }

        /// <summary>
        /// Gets the median value of a set or ordered values.
        /// </summary>
        /// <param name="array"></param>
        /// <returns></returns>
        private static double ReturnMedian(double[] array)
        {
            if (array.Length % 2 == 1)
                //if n%2 == 1 then we will have an element in the middle of the vector
                return (array[array.Length / 2]);
            else
                //otherwise we will have to calculate the average between the 2 elements between the middle position.
                //atention: the vector index starts in 0 not in 1 due to this reason it is necessary to use n/2 - 1 and n/2 instead of n/2 and n/2 +1
                return (array[array.Length / 2 - 1] + array[array.Length / 2]) / 2;
        }

        /// <summary>
        /// Calculates the max value of the data in the range of the vector.
        /// </summary>
        /// <param name="calc">
        /// Dataset
        /// </param>
        /// <returns>
        /// Max
        /// </returns>
        public double Max(double[,] calc) 
        {
            double Max = double.MinValue;

            if (calc != null)
            {
                for (int i = 0; i < calc.GetLength(0); i++)
                    for (int j = 0; j < calc.GetLength(1); j++)
                        if (!double.IsNaN(calc[i, j]) && calc[i, j] > Max)
                            Max = calc[i, j];

                return Max;
            }
            else
                return double.NaN;
        }

        /// <summary>
        /// Calculates the min value of the data in the range of the vector.
        /// </summary>
        /// <param name="calc">
        /// Dataset
        /// </param>
        /// <returns>
        /// Min
        /// </returns>
        public double Min(double[,] calc)
        {
            double Min = double.MaxValue;

            if (calc != null)
            {
                for (int i = 0; i < calc.GetLength(0); i++)
                    for (int j = 0; j < calc.GetLength(1); j++)
                        if (!double.IsNaN(calc[i, j]) && calc[i, j] < Min)
                            Min = calc[i, j];

                return Min;
            }
            else
                return double.NaN;
        }

        /// <summary>
        /// Calculates the mode value of the data in the range of the vector.
        /// </summary>
        /// <param name="calc">
        /// Dataset
        /// </param>
        /// <returns>
        /// Mode
        /// </returns>
        public double Mode(double[,] calc) 
        {
            if (calc != null)
            {
                //count the number of occorrences of a every value inside the vector.
                int[] occorrences = new int[calc.Length];
                for (int i = 0; i < calc.GetLength(0); i++)
                    for (int j = 0; j < calc.GetLength(1); j++)
                        for (int k = 0; k < calc.GetLength(0); k++)
                            for (int l = 0; l < calc.GetLength(1); l++)
                                if (calc[i, j] == calc[k, l])
                                    occorrences[i + j * calc.GetLength(0)]++;

                //returns the value that repeats itseft most.
                for (int i = 0; i < occorrences.Length; i++)
                    if (occorrences[i] == occorrences.Max() && occorrences.Max() > 1)
                        return calc[i % calc.GetLength(0), i / calc.GetLength(0)];
            }
            return double.NaN;

        }

        /// <summary>
        /// Calculates the standard deviation value of the data in the range of the vector.
        /// </summary>
        /// <param name="calc">
        /// Dataset
        /// </param>
        /// <returns>
        /// Std deviation about the mean
        /// </returns>
        public double StdDeviationMean(double[,] calc)
        {
            double std = 0;
            int n = 0;
            if (calc != null)
            {
                double mean = Mean(calc);
                for (int i = 0; i < calc.GetLength(0); i++)
                    for (int j = 0; j < calc.GetLength(1); j++)
                        if (!double.IsNaN(calc[i, j]))
                        {
                            std += Math.Pow((calc[i, j] - mean), 2);
                            n++;
                        }

                if (n > 0)
                    return Math.Sqrt(std / n);
                else
                    return double.NaN;
            }
            else 
                return double.NaN;
        }

        /// <summary>
        /// Calculates the standard deviation value of the data in the range of the vector.
        /// </summary>
        /// <param name="calc">
        /// Dataset
        /// </param>
        /// <returns>
        /// Std deviation about the mean
        /// </returns>
        public double StdDeviationMean(double[] calc)
        {
            double std = 0;
            if (calc != null)
            {
                double mean = Mean(calc);
                for (int i = 0; i < calc.Length; i++)
                    std += Math.Pow((calc[i] - mean), 2);


                return Math.Sqrt(std / (calc.Length));
            }
            else
                return double.NaN;
        }

        /// <summary>
        /// Calculates the standard deviation value of the data in the range of the vector.
        /// </summary>
        /// <param name="calc">
        /// Dataset
        /// </param>
        /// <returns>
        /// Std deviation about the median
        /// </returns>
        public double StdDeviationMedian(double[,] calc)
        {
            double std = 0;
            if (calc != null)
            {
                double median = Median(calc);
                for (int i = 0; i < calc.GetLength(0); i++)
                    for (int j = 0; j < calc.GetLength(1); j++)
                        std += Math.Pow((calc[i, j] - median), 2);


                return Math.Sqrt(std / (calc.Length));
            }
            else
                return double.NaN;
        }

        /// <summary>
        /// Calculates the coefficient of skewness value of the data in the range of the vector.
        /// </summary>
        /// <param name="calc">
        /// Dataset
        /// </param>
        /// <returns>
        /// Coefficient of skewness
        /// </returns>
        public double CoefficienfSkeness(double[,] calc) 
        {
            if (calc != null)
            {
                double std = StdDeviationMean(calc);
                double mean = Mean(calc);
                double coefficientSkeness = 0;

                for (int i = 0; i < calc.GetLength(0); i++)
                    for (int j = 0; j < calc.GetLength(1); j++)
                        coefficientSkeness += Math.Pow((calc[i, j] - mean), 3);

                return (coefficientSkeness / (calc.Length * Math.Pow(std, 3)));
            }
            else
                return double.NaN;
        }

        /// <summary>
        /// Calculates the coefficient of variation value of the data in the range of the vector.
        /// </summary>
        /// <param name="calc">
        /// Dataset
        /// </param>
        /// <returns>
        /// Coefficient of variation
        /// </returns>
        public double CoefficienfVariation(double[,] calc)
        {
            if (calc != null)
            {
                double std = StdDeviationMean(calc);
                double mean = Mean(calc);

                return std / mean;
            }
            else
                return double.NaN;
        }

        /// <summary>
        /// Calculates the interquartile range value of the data in the range of the vector.
        /// </summary>
        /// <param name="calc">
        /// Dataset
        /// </param>
        /// <returns>
        /// IQR
        /// </returns>
        public double InterquartileRange(double[,] calc) 
        {
            //To perform the interquartile range calculation I will use the median of the sub vectors
            if (calc != null)
            {
                double[] array = new double[calc.Length];

                for (int i = 0; i < calc.GetLength(0); i++)
                    for (int j = 0; j < calc.GetLength(1); j++)
                        array[i + j * calc.GetLength(0)] = calc[i, j];

                Array.Sort(array);

                if (array.Length % 2 == 0)
                {
                    double[] q1 = new double[array.Length / 2];
                    double[] q3 = new double[array.Length / 2];
                    Array.Copy(array, q1, array.Length / 2);
                    Array.Copy(array, array.Length / 2 , q3, 0, array.Length / 2);
                    return ( ReturnMedian(q3) - ReturnMedian(q1));
                }

                else 
                {
                    //if the median is a data it will be included in both halves.
                    double[] q1 = new double[array.Length / 2 + 1];
                    double[] q3 = new double[array.Length / 2 + 1];
                    Array.Copy(array, q1, array.Length / 2 + 1);
                    Array.Copy(array, array.Length / 2 + 1, q3, 0, array.Length / 2 );
                    //atention: the vector index starts in 0 not in 1 due to this reason it is necessary to use n/2 - 1 and n/2 instead of n/2 and n/2 +1
                    return (ReturnMedian(q3) - ReturnMedian(q1));
                }
            }
            else
                return double.NaN;
        }

        /// <summary>
        /// Calculates the correlation coefficient
        /// </summary>
        /// <param name="data1">
        /// First data set
        /// </param>
        /// <param name="data2">
        /// Second data set
        /// </param>
        /// <returns>
        /// Correlation coefficient
        /// </returns>
        public double CorrelationCoefficient(double[,] data1, double[,] data2) 
        {
            if (data1 != null && data2 != null && data1.Length == data2.Length)
            {
                double stddev1 = StdDeviationMean(data1);
                double stddev2 = StdDeviationMean(data2);

                return Covariance(data1, data2) / ( stddev1 * stddev2);
            }
            else
                return double.NaN;
        }

        /// <summary>
        /// Calculates the covariance
        /// </summary>
        /// <param name="data1"></param>
        /// <param name="data2"></param>
        /// <returns></returns>
        public double Covariance(double[,] data1, double[,] data2)
        {
            double p = 0;
            if (data1 != null && data2 != null && data1.Length == data2.Length)
            {
                double m1 = Mean(data1);
                double m2 = Mean(data2);

                for (int i = 0; i < data1.GetLength(0); i++)
                    for (int j = 0; j < data1.GetLength(1); j++)
                        p += (data1[i, j] - m1) * (data2[i, j] - m2);

                return p / data1.Length;
            }
            else
                return double.NaN;

        }

        /// <summary>
        /// Gets the coefficients of the line that best describes the distribution of points.
        /// </summary>
        /// <param name="data1"></param>
        /// <param name="data2"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        public void LinearRegression(double[,] data1, double[,] data2, out double a, out double b) 
        {
            // finds a, b where y = ax +b

            double xy = 0;
            double x2 = 0;
            double sumx = 0;
            double sumy = 0;


            if (data1 != null && data2 != null && data1.Length == data2.Length)
            {
                for(int i = 0; i < data1.GetLength(0); i++)
                    for (int j = 0; j < data1.GetLength(1); j++) 
                    {
                        xy += data1[i, j] * data2[i, j];
                        x2 += Math.Pow(data1[i, j], 2);
                        sumx += data1[i, j];
                        sumy += data2[i, j];
                    }

                a = (xy - (sumx * sumy) / data1.Length) / (x2 - Math.Pow(sumx, 2) / data1.Length);
                b = sumy / data1.Length - a * sumx / data1.Length;
            }
            else
            {
                a =0; 
                b = 0;
            }
        }

        /// <summary>
        /// Calculates teh correlation function of two distributions.
        /// </summary>
        /// <param name="data1"></param>
        /// <param name="data2"></param>
        /// <returns></returns>
        public double CorrelationFunction(double[,] data1, double[,] data2)
        {
            if (data1 != null && data2 != null && data1.Length == data2.Length)
            {
                double stddev1 = StdDeviationMean(data1);
                double stddev2 = StdDeviationMean(data2);

                return CovarianceFunction(data1, data2) / (stddev1 * stddev2);
            }
            else
                return double.NaN;
        }


        /// <summary>
        /// Calculates teh covariance function of two distributions
        /// </summary>
        /// <param name="data1"></param>
        /// <param name="data2"></param>
        /// <returns></returns>
        public double CovarianceFunction(double[,] data1, double[,] data2)
        {
            if (data1 != null && data2 != null && data1.Length == data2.Length)
            {
                double p = 0;

                double m1 = Mean(data1);
                double m2 = Mean(data2);

                for (int i = 0; i < data1.GetLength(0); i++)
                    for (int j = 0; j < data1.GetLength(1); j++)
                        p += (data1[i, j] * data2[i, j] - m1 * m2);


                return p / data1.Length;
            }
            else
                return double.NaN;
        }

        public double Moment_of_Inertia(double[,] data1, double[,] data2, int tolerance)
        {
            double p = 0;
            int n = 0;
            if (data1 != null && data2 != null && data1.Length == data2.Length)
            {

                for (int i = 0; i < data1.GetLength(0); i++)
                    for (int j = 0; j < data1.GetLength(1); j++)
                        if (!double.IsNaN(data1[i, j]))
                        {
                            for (int k = -tolerance / 2; k < tolerance / 2; k++)
                                for (int l = -tolerance / 2; l < tolerance / 2; l++)
                                    if ((i + k >= 0) && (j + l >= 0) && (i + k < data2.GetLength(0)) && (j + l < data2.GetLength(1)) && !double.IsNaN(data2[i + k, j + l]))
                                    {
                                            p += Math.Pow((data1[i, j] - data2[i + k, j + l]), 2);
                                            n++;
                                    }
                        }

                if (n > 0)
                    return p / (2 * n);
                else
                    return double.NaN;
            }
            else
                return double.NaN;
        }

    }
}
